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Balance equations can buffer noisy 
and sustained environmental 
perturbations of circadian clocks 
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We present a new approach to understanding how regulatory networks such as circadian 
clocks might evolve robustness to environmental fluctuations. The approach is in terms of 
new balance equations that we derive. We use it to describe how an entrained clock can 
buffer the effects of daily fluctuations in light and temperature levels. We also use it to 
study a different approach to temperature compensation where instead of considering a 
free-running clock, we study temperature buffering of the phases in a light-entrained clock, 
which we believe is a more physiological setting. 
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1. INTRODUCTION 

Circadian oscillators are entrained by the daily cycles of 
light and temperature. It is therefore important that a 
clock is sensitive to their daily periodicity. On the 
other hand, there are very substantial stochastic day- 
to-day fluctuations in these environmental cycles. 
This can, for example, be seen in the time series for 
light intensity and temperature shown in figures 1 
and 2. The daily fluctuations in both are substantial: 
the fluctuations in light have a coefficient of variation 
of approximately 36 per cent, while those of England's 
maximum and minimum temperatures in degrees centi- 
grade over just the single month of September have a 
coefficient of variation 15 and 29 per cent, respectively. 

The presence of such substantial noisy perturbations 
raises a number of questions. Our theoretical under- 
standing suggests that in a typical model these 
perturbations will produce substantial fluctuations in 
the protein levels (e.g. figure 4). This seems to be in con- 
flict with the need for the clock to provide robust signals 
to the genes that it is controlling. It therefore raises the 
question of whether the clock can be designed so that 
the daily protein variation is robustly entrained but at 
the same time the system manages to effectively 
buffer the stochastic variations. 

In this paper, we explain how it is relatively easy for 
evolution to adjust a clock network so as to carry out 
such buffering of stochastic environmental fluctuations. 
We argue that this is potentially an important reason 
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why an oscillator, rather than a direct measurer of the 
entraining signal (light or temperature) is used. The 
mathematical argument that we employ is very general 
and can be applied to a much broader class of regulat- 
ory and signalling systems and other environmental 
factors. It uses a combination of ideas behind balance 
equations (used previously to explain temperature com- 
pensation [1-3]) and the principal component aspects 
of global sensitivity analysis [4-8]. We consider fluctu- 
ations in the light level and temperature fluctuations. 
To give examples of the effectiveness of the buffering 
mechanism for light fluctuations, we use a published 
model of the Arabidopsis thaliana circadian clock [7], 
and for temperature fluctuations, we introduce a 
temperature-dependent version of this model. 

A related issue was addressed in a recent paper [8] . In 
a study involving artificial in silico evolution of clock 
networks, it was shown that a combination of the 
need to cope with multiple photoperiods and stochastic 
variation in the timing of dawn and dusk favoured the 
evolution of more loops and light inputs and greater 
complexity in the networks. 

After considering short-term fluctuations, we turn to 
temperature compensation. Temperature compensation 
refers to the striking and defining feature of circadian 
clocks whereby their period only varies by a small 
amount over a physiological range of temperatures 
[9,10]. However, the exact value of the free-running 
period in constant conditions does not appear to have 
a direct selective value in the natural environment, as 
the clock will normally be entrained to diurnal day/ 
night cycles. One may therefore ask why temperature 
compensation has arisen during evolution. We address 
this question here and develop a theory to address the 
question of how evolution might act on forced entrained 
oscillators. We show that for typical systems, the phases 
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Figure 1. Light intensity measurements for the month of 
September 2005 from the environmental radiometry data 
from Harvard Forest [17]. 



of the proteins in the clock will vary significantly with 
temperature but that one can tune the clock to satisfy 
certain balance equations so that the changes in phase 
are buffered. We propose that the buffering of the 
free-running period characteristic of classically tempera- 
ture-compensated systems is a consequence of this 
phase buffering. The idea that there is a connection 
between the free-running period and the entrained 
phase is not a new one. The way in which phase changes 
as one crosses an Arnold tongue is well understood in 
dynamical systems theory [5,11] and in circadian 
rhythms [12,13]. This relation between the period and 
the entrainment phase has been observed experimen- 
tally in physiologically relevant situations [14]. 

Throughout, we assume that our circadian clock is 
modelled by a set of ordinary differential equations 



da; 
~dt. 



= f(t, x, k), 



(1.1) 



where t is time, a vector x = (x±, . . . , x n ) represents the 
state variables (namely, mRNA and protein levels) 
and k= (ki, . . . , k s ) is a vector of parameter values. 
We assume that equation (1.1) has an attracting 
periodic solution x= g(t,k) of period T. We study 
the properties of this solution. We will illustrate our 
results by using a model of the Arabidopsis circadian 
clock [7]. 



2. GLOBAL SENSITIVITY ANALYSIS AND 
ITS PRINCIPAL COMPONENTS 

If g(t) is the solution of equation (1.1) mentioned 
above, then the change 8g(t) in g caused by a change 
8k = (Ski, . . . , 8k s ) in the parameter vector k is 




Figure 2. Daily maximum and minimum temperature 
measurements for the month of September 2005 from UK 
Met Office Historical Central England data (HadCET) [18]. 



where the linear map M is given by 
M8k = V ^-8h. 

We can regard M as a map from the parameter space W 
to a space of time series. For our purposes, the appropri- 
ate space of time series is the Hilbert space H. defined 
in appendix A, and we also consider the subspace Ho 
of H. spanned by the functions dg/dkj(t) on 0 < t < T, 
j= 1, . . . , s. 

In Rand [6] , it is shown that there are a set of num- 
bers <T\> (t 2 > . . . > <x s , a set of orthonormal vectors 
Vi, . . . ,V s oi the parameter space W and a set of ortho- 
normal vectors Ui,..., U s in H such that MVi = criUi, 
M* Ui = o~i Vi with the following optimality property: 
for all k > 1 , the average error given by 



Ml=i 



Mv ■ 



< Mv, Ui >£ 2 Ui 



dv 



is minimized over all orthonormal bases of H 0 . At this 
minimal value, ef = ccrf, where c is an absolute con- 
stant. The o"j are uniquely determined and the V t and 
Ui are, respectively, eigenvectors of MM* and M*M. 
Thus, the cr, are the eigenvalues of M*M. If they are 
simple eigenvectors, then the Ui and Vi are uniquely 
determined. 

M* is the adjoint to M and is given by 



M*U = (%,..., tj s ), 



where 



dg 



U 



It follows that the zjth element of M*Mis given by 



8g(t) = M8k+ 0(|| 8k || 2 ), 



dg dg \ 
dki ' dkjlp 
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Figure 3. Plot of log 10 (<r i /t7 , 1 ) for the largest singular values 
of the clock model [9]. The decay exponentially. 



Since this is self-adjoint, it has real positive eigenvalues 
and these are ai > cr 2 > . . . > <x 5 . 

It follows from the above discussion that a change 8k 
to the parameters leads to a change 8g in the solution g 
such that 



WijSkj 



Ui(t) + 0(|| 8k 



and that the cr^ decay as rapidly as is possible for any 
such representation. Here, the Wy are the elements of 
the inverse matrix W= V~ . The speed with which 
the o"j decay for the Arabidopsis clock model [7] is 
shown in figure 3. Inspection of figure 3 shows that 
this decay is exponential for the case for the clock 
model of Locke et al. [7], and this is typical [4,6,15] as 
also shown in the electronic supplementary material, 
figure SI. 

Note that when applying this theory, one can either 
apply it directly to the model parameters kj or one can 
apply it to the parameters r\j = log kj. It is often more 
relevant to use the second approach because, in regulat- 
ory and signalling systems, the values of two parameters 
may differ by an order of magnitude or more. Therefore, 
it is more appropriate to consider relative changes in the 
parameters kj than the absolute changes. The only 
change to the theory when considering relative changes 
instead of absolute ones is to replace the linear operator 
M above by M- A, where A = diag(fc) is the diagonal 
matrix whose diagonal is made up of the parameter 
values. If we denote quantities for the latter relative 
case using a tilde, we have V = A • V, Wy = kjW^ and 
<Ji = (Ti. This follows from the fact that for small 
changes 8k to the model's parameters, 8r)j= 8kj/kj. 
The scaled changes 8tj ; also have the advantage of 
being non-dimensional. 



3. BALANCE EQUATIONS FOR SUSTAINED 
AND DAILY FLUCTUATIONS 

Now suppose that a subset {kj t , ... , of the par- 
ameters depend upon a parameter p, i.e. fc, = kj (p). 
Then, 



dg 
dp 



where the latter sum is over j in S = {ji 
if we want dg/dp = 0, we require 



Ui{t), (3.1) 
,jg}. Thus, 



dp 



0 



(3.2) 



for i = 1, . . . , s. The s equations in (3.2) are called the 
balance equations and the sums Yljes Wy dkj/dp are 
called balance sums. 

Note that the ith balance equation can only be 
solved if, for this value of i, the Wijdkj/dp do not all 
have the same sign. This places a constraint on the 
Wy and hence on the system as a whole. 

Now we come to the reason why using the above 
principal components {/, is important. Note that in 
equation (3.1), each term J^ies Wijdkj/dp is multiplied 
by (Ti and therefore if the cr^ decrease rapidly, as usually 
is the case for the sort of systems that we consider, 
then the importance of the balance equations in ensur- 
ing dg/dp si 0 decreases rapidly as i increases, and 
obtaining the balance equations for just a few low i 
will substantially decrease |d<7/dp|. Inspection of 
figure 3 shows that this is certainly the case for 
the clock model of Locke et al. [7] and this is typical 
[4,6,15]. 

For circadian oscillators, we are often particularly 
interested in the changes in phase of the various com- 
ponents of the clock. If we define the phase <p m of the 
mth component as the time that it reaches its maximum 
value, then [6] 



d<P m 
dp 



<n{Ej<ES Wijidkj/dp)) U im (ip, 

9m(<Pm) 



(3.3) 



A derivation of equation (3.3) is provided in appendix 
A. Thus, the balance equation to obtain d(p m /dp « 0 
agrees with equation (3.2). 

In practice, the balance sums are never zero, and the 
aim of the balancing is to reduce them substantially in 
the sense that the ratio of the sum after balancing to 
that before is substantially less than one. 



4. APPLICATION TO DAILY LIGHT 
FLUCTUATIONS 

Our aim here is to demonstrate a simple mechanism 
that can enable the clock to filter out the sort of sub- 
stantial variation in daytime light intensity that is 
observed in the Harvard Forest data [17], figure 1. 

We denote the time-dependent light intensity by 
6(t). For example, 8(t) might be the function that is 1 
between dawn t\ and dusk t A and zero elsewhere or it 
might be a slightly smoothed version of this 
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Table 1. Left: parameters for the balanced model. Each light parameter kj has light intensity of the form ajL(t). For the 
Locke model, ctj= a = 1, while for the balanced model, otj= cpc + gL where Cj are listed in the table and cL are dj = 1 — Cj. 
Right: the sum Y^jWijdkj/da for i = 1,. . ., 8 evaluated for the Locke model (dfe,/da = 1) and the balanced model {dkj/da = Cj). 
The corresponding singular values cr^ are plotted in figure 3. 
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Figure 4. Effect of light intensity fluctuations on the protein 
profile of the morning loop component LHY of the Locke 
model. The profile reflects levels of the cytoplasmic protein 
for light intensity a ~ A/"(l,0.2). 

discontinuous function. We model fluctuations in the 
light by replacing 9(t) by a6(t), where a fluctuates 
daily around its mean value of a = 1. 

We then consider a model where the effect of light is 
to modulate a subset of the terms in the differential 
equation so that they are functions of a9(t). In general, 
the light input will occur in multiple terms in the form 
kj9(t), where kj is one of the parameters of the system. If 
the light level is fluctuating as above, then the term 
kj6(t) is replaced by kja6(t). Let S be the set of 
parameter indices for the parameters that occur in 
this way. 

Suppose that this original model is not buffered 
against variation in daytime light intensity. To buffer 
it, we can assume that there is some simple regulation 
of the light inputs, so each of the terms kja6(t), j G S, 
is replaced by kj(cj<x + dj)6(t). It is also natural 
to assume that c ; - + rf ; = 1 since this implies that 



kj(cja + dj) fluctuates around kj. One can always 
reduce to this case by scaling the kj in advance. We 
call this the light-modified model. 

If Cj and dj are fixed (and not regarded as par- 
ameters), then the new value for Wy in this new 
system is CjWy, where the latter Wy is the value in 
the original system. Therefore, the balance equation 
for absolute changes in the parameters is 

o", W v c j k J = °~> ^« °i = °> ( 4 " 1 ) 

where the are as above for the approach where one 
uses relative changes in parameters instead of absolute 
ones so that one uses rjj = log kj. 

This is a very general formulation and is, for example, 
directly applicable to the model of Locke et al. [7] for the 
Arabidopsis clock that introduces light in the way 
described. The light parameters kj, j G S, are listed in 
table 1. The light parameters that have the greatest 
effect on the balance equations (by ranking o-iWy) 
are LHY transcription (q^, light and TOC1- 
mediated induction (g 2 and n 4 ) of Y transcription and 
accumulation of protein P (p 5 ). 

We assume that the variations a are normally dis- 
tributed with mean /jl and standard deviation cr, 
which we denote by a ~A/"(/a,o"). As shown in figure 4 
(and in the electronic supplementary material, figure 
S2), the daily variation in the amplitude of the limit 
cycle solution of this model is substantial and quanti- 
tatively reflects the variation in the light amplitude. 
A light-modified model was constructed as described 
for which the left-hand sides of equation (4.1) are 
substantially smaller, as is given in table 1. 

To implement this, the values of six of the Cj were 
chosen fairly arbitrarily but relatively close to one and 
then the remaining three were calculated by solving 
the linear balance equations corresponding to the first 
three principal components. Our balanced model 
shows less variation in the output protein and gene 
levels (figure 5 for a ~ Af(l,0.2)). For all genes and 
proteins, the balanced model shows less variation in 
the concentration levels than the unbalanced Locke 
model (electronic supplementary material, figure S3 
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Figure 5. Effect of light intensity fluctuations on the morning 
loop protein LHY of the balanced model. The profile reflects 
levels of the cytoplasmic protein for light intensity a ~ 
Af(l,0.2). 



and table SI). The phases of the balanced model show 
less variation for all components, except the LHY 
gene and proteins, which are of same order as those of 
the Locke model (electronic supplementary material, 
table S2). 



5. APPLICATION TO DAILY 

TEMPERATURE FLUCTUATIONS 

We consider a model with fixed but differing day and 
night temperatures, T D and T N . Below, we derive a 
set of balance equations that eliminate the effects of 
temperature fluctuations. 

We assume that the temperature enters the model 
through a set of parameters kj, j £ S, which are temp- 
erature sensitive. A standard assumption for the 
temperature dependence of each of model parameters 
kj is that it is similar to that for rate constants of chemi- 
cal reactions and is described by the Arrhenius 
equation. This expresses the dependence of the rate con- 
stant kj on the temperature T and activation energy Ej 
as kj= Ajexp(— Ej/RT), where A, is a constant specific 
to the individual parameter and R is the gas constant 
(8.314472 x 10" 3 kj mol^K" 1 ). In our case, we need 
to deal with the fact that we have different night and 
day temperatures and therefore we assume 



kj(t) 



Aj exp 



Aj exp 



■EL 



R(T D + e) 
R(T N + V ) 



if t is in the daytime 

if t is in the night, 

(5.1) 



where t is in the daytime if t\< t mod r < t&, where 
[t h £ d ] denotes the range of day hours, and otherwise 
t is in the night. Moreover, e and r\ denote the 



fluctuations in day and night temperatures, T D and 
T N , respectively. We consider the following first-order 
Taylor series expansion: 

&7 n + fc?, e if £ is in the daytime 
kj 0 + k } - x 17 if t is in the night, 

where fc° 0 = Aj exp ( — Ej/ R T D ) and kf A = - k^ 0 E } /RT u 
and similarly for the night parameters. Note that this is 
a very good approximation and higher order terms can 
be neglected, since for a parameter of order O(0. 1) (i.e. 
the order of a large number of Locke parameters), with 
sensible activation energy 50 kJ mol" 1 , the coeffi- 
cient of the second-order term is of order O(0. 001). 

From the observation in equation (3.1), day temp- 
erature variations e will have the following effects on 
changes to the solution g: 



dg 
de 



E°nE m^iW) 



(5.3) 



and the following changes to the phases: 



d<P„ 



de 



^{Zjes " / -':J UiteJ 

9 m (<Pm) 



(5.4) 



Similar expressions can be derived for night tempera- 
ture fluctuations. Together, these give rise to two sets of 
balance equations, 



<^E 



0 and 



0. (5.5) 



Or alternatively, 

T^E W v^ = 0 and -^Yl 

D jfES N jes 



W^Ej = 0, 



(5.6) 



since, by the above, K?i = — kfoEj/ ' RT\\. In order to 
make a model temperature dependent, we have to 
define the dependence of its parameters on temperature. 
It is likely that all parameters in a regulatory system 
are actually temperature dependent, but this tem- 
perature dependence will have little effect for a 



parameter kj if all the sensitivities 5L- : 



<x, : Wij are 



small for that value of j. Therefore, in order to keep 
the model reasonably tractable, we will only introduce 
temperature into those variables with a significant 
sensitivity. 

Moreover, to determine the relative importance of 
parameters, since values of some parameters differ by 
an order of magnitude or more, it is more appropriate 
to compare relative changes of parameters and to use 
the log parameters and the Wy as described in the 
last paragraph of §2. We selected parameters that 
ranked highest when ordered by max i=1 _ 4 |o-jW^|. 
The plot of top 20 is shown in figure 6. In this selection, 
we include only parameters that come linearly in the 
model and exclude Hill coefficients as these are not 
rate constants. Parameter has the highest max i=1 _ 4 
logio| cTjlly. Only seven other parameters have 
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Figure 6. Power sensitivity spectrum of the Locke model. Each group of bars corresponds to the values of log 10 |(TjW^| for a 
parameter kj. These are only plotted for those i for which logio|(7i is significant (here that is i= 1-4). The parameters kj 
are ordered by max j=1 _ 4 logiolcTjWJ and only 20 most sensitive ones are plotted. Not taking into account Hill terms g 2 and 
g-j (for reasons outlined in the main text), parameter p e has the highest max i=1 _ 4 logrol (TjWy. Only seven other parameters (par- 
ameters up to and including mi) have sensitivity higher than 30 per cent of the maximum, and they will be candidates for 
temperature-sensitive parameters. 

Table 2. Energies and night and day parameter values for the unbalanced and the balanced model. Note that day temperature 
r D =290.15 K and night temperature is T N =285.15 K. 



unbalanced balanced 



kj 




$J (To) 


(Tn) 


E 3 




koj(Tn) 


r>e 


13.7721 


7.6670 


8.4742 


27.6101 


7.2635 


8.8777 


m m 


39.1133 


11.0158 


14.6380 


27.6091 


11.0158 


13.4638 


Pe 


41.5837 


0.2471 


0.3343 


27.6374 


0.2616 


0.3198 


m 17 


2.7490 


4.4595 


4.5495 


2.7490 


4.4595 


4.5495 


m 4 


13.7658 


3.6320 


4.0142 


27.6079 


3.4408 


4.2054 


n 0 


27.5836 


0.1991 


0.2433 


27.5836 


0.1991 


0.2433 




27.6115 


2.7078 


3.3096 


27.6115 


2.7078 


3.3096 


mi 


27.6087 


1.7991 


2.1989 


27.6087 


1.7991 


2.1989 





pc 1 




pc 2 




pc 3 




pc 4 



sensitivity higher than 30 per cent of the maximum. 
These eight parameters will be temperature sensitive. 

Note that the key parameters to control the buffering 
were parameters linked to PRR7/9 components, 
namely LHY-dependent transcription (rzg), mRNA 
degradation (m 16 ), translation (p@) and cytoplasmic 
protein degradation (m 17 ) as well as the nuclear- 
cytoplasmic transport (ri 2 ). Aside from these, the list 
also includes TOC1 light-independent transcription 
(n 2 ), mRNA degradation (77J4) and LHY mRNA degra- 
dation (mi). Each of the parameters kj was separated 
into a morning and an evening term, kj=kf9 + 
kf{l — 9), where morning and evening terms were 
adjusted to be a fixed percentage higher or lower than 
the Locke model values. Note that the term 6 is the 
light function. For most parameters, we chose this to 
be 10 per cent, but for 77117, we had to choose a signifi- 
cantly smaller value of 1 per cent, so that gene and 



protein waveforms of the temperature-sensitive model 
would match those of the Locke model as closely as 
possible. We confirmed that the parameters had a 
strong effect on the model by verifying that several of 
them appeared in the top 10 per cent of parameters 
ordered by their magnitude of maxj =1 _4logi 0 |o"j Wd. 
We chose night and day temperatures to be T = 
285.15 K and T= 290.15 K, to be close to the tempera- 
ture data from UK Met Office (figure 2), with mean 
maximum and minimum temperatures T= 292.79 K 
and T= 283.92 K. 

We calculated the energies of the new model and 
checked the balance equations. We tried several differ- 
ent combinations of energies, then recalculated the 
morning and evening parameters and then computed 
the balance equations for the model. From the starting 
Locke model, we evolved two models with differing 
values of energies (table 2). We labelled them a 
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Table 3. The balance sums YljWijKoj Ej an d Xw^HW ano - the corresponding singular values cr^ for the unbalanced model 
(UB) and the balanced model (B). To get the true sums, divide each column by 1/RT 2 , where r D =290.15K and 
T N =285.15 K. 
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Figure 7. Effect of night temperature fluctuations on protein 
profile of the morning loop component LHY of the Locke 
model. The profile reflects levels of the cytoplasmic protein 
for night temperature variations 17, e ~ JV(Q,1). 

balanced and an unbalanced model according to their 
fit to the balance equations (table 3). 

The balanced model was chosen so that originally it 
had activation energies of about 30 kJ mol -1 for almost 
every component except for that of the parameter m 17 
whose activation energy was chosen to be significantly 
lower to make balancing easier. Since later we adjusted 
the temperature range to fit with the data seen in 
figure 2, these values appeared slightly smaller. The 
unbalanced model was made by changing four energies 
from the list of the balanced model, in order to make a 
worse fit to the balance equations. 

In fact, our initial temperature-sensitive model gave 
the best balance equations, so it was chosen as the 
balanced model. The balanced model shows less vari- 
ation in peak concentrations and phase variations 
than the unbalanced model (figures 7 and 8 and elec- 
tronic supplementary material, tables S3 and S4). 



6. APPLICATION TO TEMPERATURE 
COMPENSATION 

Instead of temperature compensation on a free-running 
clock (i.e. dp/dTK, 0), we are interested in temperature 
compensation of the light-entrained clock, which we 
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Figure 8. Effect of night temperature fluctuations on the 
morning loop protein LHY of the balanced model. The profile 
reflects levels of the cytoplasmic protein for night temperature 
variations 17, e ~ 7V(0,1). 

believe is a more physiological setting. The aim is to 
minimize protein and phase changes in the context 
of changing temperature, dtp/dTrz 0. We therefore 
consider models that are entrained to the day-night 
cycle rather than the free-running clock. 

This hypothesis [16] that a balance of opposing reac- 
tions could allow temperature compensation was first 
tested by Ruoff using a simple model for an oscillatory 
feedback loop [1-3]. He used an Arrhenius represen- 
tation for temperature dependence and deduced a 
balance equation for the local period slope dp/dT in 
terms of the activation energies Ej and control 
coefficients for each of the parameters. 

We also assume that parameters fc, j G S, are temp- 
erature dependent and describe them by Arrhenius 
equations, kj= Aj exp(—Ej/RT), as described above. 
The temperature t is in kelvins in the range 285.15 
K < T< 300.15 K. Activation energies, Ej, must be in 
the range 1 kJ mol -1 < Ej< 150 kJ mol -1 . We insist 
that some of the activation energies are substantial 
because otherwise the parameters only have weak depen- 
dence upon temperature. In fact, for a given balanced 
system, since these energies enter linearly into the bal- 
ance equation, scaling them by a factor just scales the 
divergence from perfect balance by that factor. 
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Table 4. Energy values Ej of the temperature-sensitive 
models based on the Locke et al. [7] model: model 1 (Ml) 
and model 2 (M2). 



kj 


Ej (Ml) 


Ej (M2) 


n& 


11.3399 


22 


m m 


12.3650 


5 


Pa 


7.5775 


7.5775 




10.2539 


10.2539 


no 


48.3543 


48.3543 


m 4 


42.2986 


10 


fh 


50.8636 


50.8636 


m\ 


1.2201 


1.2201 



We now apply equation (3.2) where p is replaced by 
temperature T. From the relation dkj /dT = kjEj / RT 2 1 
we deduce that the balance equations are 

for i= 1, . . ., s. When the cr, decrease rapidly, we need 
only consider these balance equations for the first few 
i in order to get dg/dT or d<p m /dT small. Note that 
the ith equation can only be solved if, for this value of 
i, the Wy do not all have the same sign. This is because 
kj Ej is always positive. This places a constraint on the 
and hence on the system as a whole. Only certain 
networks can be balanced. 

We assume that the model parameters from Locke 
et al. [7] correspond to a model at T= 288.15 K. Since 
this model does not include temperature, we selected 
the temperature-sensitive parameters at each tempera- 
ture end by checking the parameter sensitivity 
spectrum, as described in §5. 

We consider three models in which temperature 
dependence is inserted into the Locke model. In two 
of them (models 1 and 2), the parameter values at 
T= 288.15 K are as in the original model but although 
they have the same activation coefficients Aj, they have 
different activation energies Ej (table 4). 

The energies for model 1 were selected so that they 
were substantial and roughly satisfied the balance 
equation (table 6). We then iteratively calculated the 
Wijs for the new model and more exactly rebalanced 
the equation. The iterative procedure converged reason- 
ably well. For model 2, the energies were modified so 
that the balance equations were not well satisfied 
(table 6). Some energies were decreased as well as 
increased (table 4). We could easily do this as the two 
models share W^s at T= 288.15 K. 

The temperature compensation in model 2 is sub- 
stantially worse than that in model 1, demonstrating 
the importance of balancing. We find that model 1 
is better at local temperature compensation than 
model 2, with less variation in its protein concentrations 
and phases from the original published model [7], 
cf. TOC1 protein times series in figure 9, and also 
it fares better at global temperature compensation 
(electronic supplementary material, tables S8 and S9). 
Moreover, model 1 can temperature compensate in con- 
stant light conditions (figure 10), with model 2 
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Figure 9. TOC1 protein profiles. The Locke model at T = 
288.15 K (dashed line) against the balanced model 1 at T = 
289.15 K (a) and the unbalanced model 2 at T= 289.15 K 
(b). By balancing, we can ensure minimum change to protein 
profile and phases. 

becoming arrhythmic outside a narrow temperature 
range. Model 3 shares activation energies of model 1 
but has different activation coefficients Aj. Thus, its 
parameter values do not agree exactly with the Locke 
model at any temperature (table 5). This model is 
better balanced than either of the other two (table 6) 
and is significantly better compensated for phase 
than the other models (tables 7 and 8 and electronic 
supplementary material, table S10). Moreover, in 
continuous light conditions, its period is much better 
compensated than model 2 and slightly better than 
model 1 (figure 10). 



7. DISCUSSION 

We have attempted to show how a combination of ideas 
behind balance equations and the principal component 
aspects of global sensitivity analysis gives a new 



Interface Focus (2011) 



Balance equations and circadian clocks M. Domijan & D. A. Rand 185 




temperature (K) 



Figure 10. Period of the free-running temperature-sensitive 
models. The plot shows the period of model 3 (open circles) 
and the Locke temperature-sensitive model with different ener- 
gies, models 1 and 2 (black squares and open triangles, 
respectively). The original model of Locke et al. [7] is plotted 
as the third black square (T= 288.15 K). Both temperature- 
sensitive Locke models cannot temperature compensate as well 
as model 3. Model 2 is arrhythmic outside the range plotted. 



Table 5. Values of temperature-dependent parameters for 
model 3 with T 0 = 285.15 K and 2\ = 300.15 K. 



kj 


E 3 


h (To) 


^ (2\) 




11.3399 


7.6784 


9.7517 


m 16 


12.3650 


11.5927 


15.0445 


Pe, 


7.5775 


0.2812 


0.3299 


m 17 


10.2539 


4.2545 


5.2810 


rio 


48.3543 


0.2050 


0.5000 


m 4 


42.2986 


3.0918 


8.5677 


"2 


50.8636 


2.4062 


7.0300 


OTi 


1.2201 


1.9883 


2.0401 



approach to understanding how to design regulatory 
networks that are buffered against certain fluctuating 
environmental perturbations. As well as presenting 
some concrete examples of applications to circadian 
clocks, we have described the general theory behind 
this. From this, it is clear that it can be applied to a 
much broader class of regulatory and signalling systems 
and environmental factors other than light and 
temperature. 

If buffering of a particular environmental fluctuation 
has significant selective pressure for the organism in 
question, then we can interpret this selective pressure 
as acting on the balance equation. There will then be 
selective pressure on the quantities that make up the 
left-hand side of the equation (e.g. system parameters, 
activation energies) to make them balance to zero. 
Understanding this makes it much clearer how evolution 
can act to achieve what appear to be quite complex tasks. 

Temperature compensation has been one of the driving 
dogmas of circadian biology and has been interpreted in 
terms of the constancy or near-constancy of the free- 
running period of the circadian clock under changing 
temperature. However, it is not clear how evolution acts 



Table 6. Sums X^eS^sA-^j at each temperature T a = 
285.15 K, T x = 288.15 K and T 2 = 300.15 K for all three 
models. To get the true sums, divide each sum by 1/RT 2 . The 
corresponding singular values a are shown in parentheses. 



J2j eS W&fflEj and fa (x 10 4 )) 





i T 0 T x 


T; 


model 1 


1 0.0397 


0.0419 


-0.0229 




(18.8826) 


(1.0648) 


(1.2884) 




2 -1.1872 


0.0134 


0.1168 




(0.0261) 


(0.0352) 


(0.0367) 


model 2 


1 0.0840 


-0.0787 


0.0544 




(4.8274) 


(1.0648) 


(3.6128) 




2 1.0042 


-0.2398 


0.1579 




(0.0325) 


(0.0352) 


(0.0310) 


model 3 


1 0.0434 


0.0286 


0.0162 




(1.2388) 


(1.5245) 


(1.2037) 




2 -0.0681 


0.1220 


0.1057 




(0.0368) 


(0.0351) 


(0.0368) 


Table 7. 


Peak and trough times of 


cytoplasmic protein at 


temperatures T a = 285.15 K and T\ 


= 300.15 K of model 3 


(the globally compensated model) 


with parameter values 


from table 5. 
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(T= T 0 ) 
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x c 


5.0 4.5 


15.7 


15.7 


Y c 
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PRR7/9, 


16.9 16.9 
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Table 8. 


Peak and trough times of 


cytoplasmic protein of 


model 2 at temperatures T a = 285.15 K and 2\ = 


300.15 K. 




peak times 


trough times 




(T= T 0 ) (T= T%) 


(T= T 0 ) 


(T= J\) 


LHY C 


7.1 8.2 


19.8 


21.9 


TOCl c 


23.1 13.1 


9.5 


23.5 


x c 


5.3 20.1 


15.9 


5.0 


Y c 


6.9, 16.4 2.2, 7.0 


2.4, 9.8 


6.0, 17.5 


PRR7/9, 


16.0 16.9 


0 


2.5 



on the free-running period since in physiological con- 
ditions the clock is entrained to the day-night cycle 
and has a period of 24 h. We show that through certain 
balance equations it is possible to buffer the changes in 
phase over relatively large ranges of temperature. 
Although we do this for a specific example, the math- 
ematical approach suggests that this buffering should 
be possible for an extremely broad range of clock 
models. We suggest that the observed near-constancy 
of the free-running period is a consequence of the near- 
constancy of the phases of the entrained clock. 

To balance a balance equation, it is necessary that 
the terms making up the equation do not all have the 
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same sign. For example, for temperature compensation, 
we require that the relevant quantities Wy do not all 
have the same sign. This puts constraints on the net- 
work structure, and this leads to a prediction about 
what network structures can be expected. 



8. METHODS 

All the computations were carried out using Matlab and 
XPPAUT. In particular, the global sensitivity calcu- 
lations were done using the Matlab-based Time Series 
Sensitivity Analysis Package available from http:// 
www2.warwick.ac.uk/fac/sci/systemsbiology/software/, 
and the period calculations were performed using 
XPPAUT, available from http://www.math.pitt.edu/ 
~bard /xpp/xpp . html. 

We are particularly grateful to Paul Brown who contributed 
extensively to the software used by us in preparing this 
paper. We are also grateful to the ROBuST team for many 
useful discussions on this topic and particularly to Andrew 
Millar. The comments of two anonymous referees were very 
helpful. This research was funded by BBSRC SABR grant 
BB/F005261/1 (ROBuST Project) and EU BIOSIM 
Network Contract 005137. D.A.R. is also funded by EPSRC 
Senior Research Fellowship EP/C544587/1. 



APPENDIX A 

A. 1. Definition of Hilbert space H 

This is the L? Hilbert space of Revalued functions 
U(t) = (U^t), . v ,u n (t)), U'(t) = (if^t), . . .,U n (t)), 
0 < t < T, with inner product 



{U,U') P = T- 



and norm given by 



T n 



A. 2. Derivation of equation (3.3) 

This was presented in Rand [6]. Let t= <p m (k) be the 
time when the concentration of the mth component, 
g m , is at a maximum or a minimum value. Hence, 
4>m{k) satisfies g m (4> m (k),k) = 0. Differentiating both 
sides with respect to kj and rewriting, 

g^m = {dg m /dkj)(<j> m (k),k) 

9kj g,n(<t>m) 

(d/dt^JdgJdk^Jk)^) 
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